rm(list = ls())

suppressMessages(library(tidyverse))
suppressMessages(library(magrittr))

c(
list.files("/Pub/Users/wangyk/Project_wangyk/Codelib_YK/some_scr", full.names = T, recursive = T, pattern = "\\.R$"),
list.files("/Pub/Users/cuiye/RCodes/UserCode", recursive = T, full.names = T, pattern = "\\.R$")) %>% 
walk(source)

library(Seurat)
source("/Pub/Users/wangyk/Project_wangyk/Codelib_YK/some_cancers/major_test2.r")
# conflicted::conflict_prefer_all("dplyr", quiet = T)

out_home <- "/Pub/Users/wangyk/project/Poroject/F230731001_OV_sc/"
setwd(out_home)

t_sub <- readRDS("out/SC_part/04.tsub2/t_sub_seuratObj.rds")

library(CellChat)
library(patchwork)
options(stringsAsFactors = FALSE)

# t_sub@assays$RNA@data
t_sub@meta.data %>% select(CellStatus) %>% head()
table(t_sub@meta.data$BRCA_mut)

# 做三次，分组以及整体---------
m1 <- subset(t_sub,BRCA_mut == "BRCA1_mut")
m2 <- subset(t_sub,BRCA_mut == "BRCA2_mut")
m12 <- subset(t_sub,BRCA_mut != "WT")

# https://rpubs.com/DarrenVan/879013 非常好的教程

cellchat_simple <- \(seurat_obj,col_in_meta = 'CellStatus',od,var_name){
    # 构建对象
    cellchat_obj <- createCellChat(
        object = seurat_obj@assays$RNA@data,
        meta = seurat_obj@meta.data %>% select(any_of(col_in_meta)),
        group.by = col_in_meta
    )
    cellchat_obj@DB <- CellChatDB.human

    cellchat_obj <- setIdent(cellchat_obj, ident.use = col_in_meta)
    cellchat_obj <- subsetData(cellchat_obj)

    # subset the expression data of signaling genes for saving computation cost

    # 设置并行计算，设置每个worker的最大内存
    max_memory <- 4000
    options(future.globals.maxSize = max_memory * 1024^2)
    future::plan("multiprocess", workers = 20)
    # 计算过表达的signaling genes and LR pairs
    message("Calculating for over-expressed signaling genes associated with each cell group ..." %>% cli::col_blue())
    cellchat_obj <- identifyOverExpressedGenes(cellchat_obj)
    message("Calculating for over-expressed ligand-receptor interactions (pairs) within the used CellChatDB ..." %>% 
    cli::col_blue())
    cellchat_obj <- identifyOverExpressedInteractions(cellchat_obj)
    raw.use <- TRUE
    population.size <- FALSE
    # 推断细胞间相互作用
    message("Working in computeCommunProb..." %>% 
    cli::col_blue())

    cellchat_obj <- computeCommunProb(
        object = cellchat_obj, type = "triMean",
        population.size = F,
        raw.use = T, do.fast = TRUE
    )

    # computeCommunProb将返回CellChat对象，该对象将包含slot 'net'
    # object@net$pval is the corresponding p-values of each interaction
    # 过滤只有少数细胞的互作结果
    cellchat_obj <- filterCommunication(cellchat_obj, min.cells = 10)
    # Compute the communication probability on signaling pathway level by summarizing all related ligands/receptors
    cellchat_obj <- computeCommunProbPathway(cellchat_obj, thresh = 0.05)
    cellchat_obj <- aggregateNet(cellchat_obj)
    # Compute the network centrality scores
    cellchat_obj <- netAnalysis_computeCentrality(cellchat_obj, slot.name = "netP")
    # 退出并行
    future::plan("sequential")

    dirname(str_glue("{od}/{var_name}/cellchat_obj.rds")) %>% mkdir

    saveRDS(cellchat_obj, str_glue("{od}/{var_name}/cellchat_obj.rds"))

    groupSize <- as.numeric(table(cellchat_obj@idents))
    pdf(str_glue("{od}/{var_name}/cell_chat_count.pdf"), width = 6, height = 6, onefile = T)
    netVisual_circle(cellchat_obj@net$count, vertex.weight = groupSize, weight.scale = T, label.edge = F, title.name = "Number of interactions")
    dev.off()

    pdf(str_glue("{od}/{var_name}/cell_chat_weight.pdf"), width = 6, height = 6, onefile = T)
    netVisual_circle(cellchat_obj@net$weight, vertex.weight = groupSize, weight.scale = T, label.edge = F, title.name = "Interaction weights/strength")
    dev.off()
}

# walk2(
#     c(m1, m2, m12), c("BRCA1_mut", "BRCA2_mut", "both"),
#     ~ cellchat_simple(seurat_obj = .x, od = "out/SC_part/05.T_sub_cellchat", var_name = .y)
# )

walk2(
    c(m2), c( "BRCA2_mut"),
    ~ cellchat_simple(seurat_obj = .x, od = "out/SC_part/05.T_sub_cellchat", var_name = .y)
)

ss_seuratObj <- readRDS("out/SC_part/02.anno/ss_seuratObj.rds")

ss_seuratObj@meta.data %<>% mutate(BRCA_mut = case_when(
    orig.ident == "GSM5599225_Cancer1" ~ "BRCA2_mut",
    orig.ident == "GSM5599229_Cancer5" ~ "BRCA1_mut",
    T ~ "WT",
))

daqun_m1 <- subset(ss_seuratObj,BRCA_mut == "BRCA1_mut")
daqun_m2 <- subset(ss_seuratObj,BRCA_mut == "BRCA2_mut")
daqun_m12 <- subset(ss_seuratObj,BRCA_mut != "WT")


walk2(
    c(daqun_m1, daqun_m2, daqun_m12), c("all_cell_BRCA1_mut", "all_cell_BRCA2_mut", "all_cell_both"),
    ~ cellchat_simple(seurat_obj = .x,col_in_meta = "CellType", od = "out/SC_part/05.T_sub_cellchat", var_name = .y)
)

list.dirs("out/SC_part/05.T_sub_cellchat")[-1]

obj_files_dir <- str_c(list.dirs("out/SC_part/05.T_sub_cellchat")[-1],"/cellchat_obj.rds")

f <- \(x){
    x = obj_files_dir[3]
    cell_chara <- ifelse(str_detect(x,"all_cell"),"T cell","CD8+ Trm")

    cellchat_obj_res <<- readRDS(x)
    df.net <- subsetCommunication(cellchat_obj_res)
    # pdf(paste0(dirname(x), "/bubble_cellchat.pdf"), width = 4.5, height = 9)
    a <- netVisual_bubble(cellchat_obj_res,
        sources.use = cell_chara, targets.use = unique(df.net$source),
        remove.isolate = FALSE, color.heatmap = c("Spectral"),return.data = T
    )

    # pdf(paste0(dirname(x), "/bubble_cellchat2.pdf"), width = 4.5, height = 9)
    a2 <- netVisual_bubble(cellchat_obj_res,
        sources.use =unique(df.net$source), targets.use =  cell_chara,
        remove.isolate = FALSE, color.heatmap = c("Spectral"),return.data = T
    )

    p <- ggpubr::ggarrange(a$gg.obj,a2$gg.obj,nrow = 1,ncol = 2,common.legend = F)
    ggsave(plot = p,filename = paste0(dirname(x), "/bubble_cellchat.pdf"), width = 9, height = 9)
}

walk(obj_files_dir,f)



# GSEA
m12 <- subset(t_sub,BRCA_mut != "WT")

c1 <- read.delim("out/SC_part/05.T_sub_cellchat/gsea/find_all_marker.xls")$cluster %>% unique() %>% .[7]
c2 <- read.delim("out/SC_part/05.T_sub_cellchat/gsea/find_all_marker.xls")$cluster %>% unique() %>% .[-7]

res <- FoldChange(m12, ident.1 = c1, ident.2 = c2)

deg_all <- res %>%
    arrange(desc(avg_log2FC)) %>%
    dplyr::select(avg_log2FC) %>% 
    rownames_to_column('gene') %>% 
    dplyr::rename(log2FC = 2)

write_tsv(deg_all,file = "out/SC_part/05.T_sub_cellchat/gsea/deg_all.txt")

# /Pub/Users/fuxj/.conda/envs/R_clusterProfiler/bin/R  
suppressMessages(library(tidyverse))
suppressMessages(library(magrittr))
deg_all <- read.delim("out/SC_part/05.T_sub_cellchat/gsea/deg_all.txt")
gene_name <- clusterProfiler::bitr(deg_all$gene,
    fromType = "SYMBOL", toType = "ENTREZID",
    OrgDb = "org.Hs.eg.db", drop = F
)

tmp_data <- left_join(gene_name, deg_all %>% dplyr::rename(SYMBOL = gene))
tmp_data <- na.omit(tmp_data)

genelist <- tmp_data$log2FC
names(genelist) <- tmp_data$ENTREZID

gsea_kegg <- clusterProfiler::gseKEGG(
    geneList = genelist,
    verbose = T,
    seed = F,
    # nPerm = 1000,#
    keyType = "kegg", # 可以选择"kegg",'ncbi-geneid', 'ncib-proteinid' and 'uniprot'
    organism = "hsa", # 定义物种,
    pvalueCutoff = .05, # 自定义pvalue的范围
    pAdjustMethod = 'BH' # 校正p值的方法
)

saveRDS(gsea_kegg,file = "out/SC_part/05.T_sub_cellchat/gsea/gsea_kegg.rds")


gsea_kegg <- readRDS("out/SC_part/05.T_sub_cellchat/gsea/gsea_kegg.rds")
gsea_kegg %>% as.data.frame() %>% .[["Description"]]

write_tsv(gsea_kegg %>% as.data.frame() ,file = "out/SC_part/05.T_sub_cellchat/gsea/gsea_kegg.xls")

f <- \(x){
    p_gsea_kegg_line <- enrichplot::gseaplot2(
        x = gsea_kegg,
        geneSetID = gsea_kegg$ID[x], # 只显示前4个GSEA的结果
        title = "KEGG", # 标题
        color = "#4DAF4A", # 颜色
        pvalue_table = FALSE,
        ES_geom = "line",
        rel_heights = c(1, 0.2, .5)
    )
    df <- as.data.frame(gsea_kegg[x, ])
    p_chara <- ifelse(df$p.adjust < 0.001, "< 0.001", sprintf(" = %.3f", df$p.adjust))

    y_pos <- ifelse(df$NES > 0,.2,-.1)

    annote_chara <- str_glue("{df$Description}\nNES = {sprintf(\"%.3f\",df$NES)}\np.adjust{p_chara}")
    p_gsea_kegg_line[[1]] <- p_gsea_kegg_line[[1]] + annotate(
        geom = "text", x = 500, y = y_pos,
        label = annote_chara, size = 4, hjust = "left"
    )

    plotout(
        p = p_gsea_kegg_line, od = paste0("out/SC_part/05.T_sub_cellchat/gsea"),
        name = str_glue("gsea_{gsea_kegg$Description[x]}"), w = 5, h = 4.2
    )

    return(p_gsea_kegg_line)
}

p_lst <- map(c(25,23), f)

p <- ggpubr::ggarrange(plotlist = p_lst, nrow = 2, ncol = 1)
plotout(
    p = p, od = paste0("out/SC_part/05.T_sub_cellchat/gsea"),
    name = str_glue("gsea"), w = 4.2, h = 8
)


